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Unbiased stochastic sampling of the one- and two-body reduced density matrices is achieved in 
full configuration interaction quantum Monte Carlo with the introduction of a second, “replica” 
ensemble of walkers, whose population evolves in imaginary time independently from the first, and 
which entails only modest additional computational overheads. The matrices obtained from this 
approach are shown to be representative of full configuration-interaction quality, and hence provide 
a realistic opportunity to achieve high-quality results for a range of properties whose operators 
do not necessarily commute with the hamiltonian. A density-matrix formulated quasi-variational 
energy estimator having been already proposed and investigated, the present work extends the 
scope of the theory to take in studies of analytic nuclear forces, molecular dipole moments and 
polarisabilities, with extensive comparison to exact results where possible. These new results confirm 
the suitability of the sampling technique and, where sufficiently large basis sets are available, achieve 
close agreement with experimental values, expanding the scope of the method to new areas of 
investigation. 


INTRODUCTION 


The full configuration interaction quantum Monte 
Carlo method (FCIQMC) and its initiator adaptation 
(j-FCIQMC) are projector QMC techniques, capable of 
providing near-exact, systematically improvable descrip¬ 
tions of correlated wavefunctions expressed as linear com¬ 
binations of Slater determinants. [l|, [2| This convergence 
is achieved by stochastically sampling the exponentially 
large (though finite) Hilbert spaces of configuration inter¬ 
action theory via a population dynamics performed on an 
ensemble of signed walkers. Annihilation processes pro¬ 
vide a means of combating the ill effects of the fermion 
sign problem which plagues projector approaches, 
exploiting the sparsity of the wavefunction induced by 
a judicious choice of orbital basis. The approach re¬ 
quires substantially less computational effort than an 
iterative diagonalisation technique, and has thus found 
considerable success in studies of atomic and molecular 
systems, [iliH model systems such as the homogeneous 
electron gas and the Hubbard model, [iMi and solid- 
state systems, [l^ 

The principal focus of many of these studies has 
been to derive properties based upon total energies, 
for which an unbiased projected estimator is readily 
available, and which have included excitation and dis¬ 
sociation energies, electron affinities, Q ionisation 

potentials, [ 1 , 121 and equations of state, [l^ Despite their 
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success, however, the extension to include the calcula¬ 
tion of a greater range of properties — expectation val¬ 
ues of operators which do not necessarily commute with 
the hamiltonian — remains highly desirable. This focus 
has been the subject of considerable interest for the QMC 
community in general, and has posed a considerable chal¬ 
lenge for decades. [l7H29| These properties, which include 
static correlation functions and entropy estimators as 
well as the forces, multipole moments, and polarisabil¬ 
ities considered here, may be deduced from the effect of 
a perturbation from the corresponding operator, P, upon 
the hamiltonian, 

H' = H + XP, (1) 


with A the perturbation strength, such that the expecta¬ 
tion value, (P), is given by the derivative of the energy 
with respect to A, evaluated at A = 0: 


{P) 


dE {H') 


dX 


A=0 


( 2 ) 


In accordance with the Hellmann-Feynman theorem, 
applicable to converged (normalised) z-FCIQMC wave- 
functions by analogy with deterministic and strictly vari¬ 
ational FCI, this expression reduces to 


(P) = {A>\P\A>), (3) 


or equivalently to the trace of P with the appropriate 
rank of reduced density matrix. ^ It is worth noting 
that unconverged f-FCIQMC wavefunctions need not rig¬ 
orously obey the Hellmann-Feynman theorem, and so in 
this work we ensure that we are working in the large 
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walker limit, such that systematic errors in the sampled 
distribution due to insufficient walker numbers have been 
minimized to the FCI-limit. 

The effective stochastic acquisition of these reduced 
density matrices, therefore, has the capacity to broaden 
the scope of i-FCIQMC significantly, and motivates the 
present work. We begin with a brief overview of the 
j-FCIQMC algorithm, including its extension to non¬ 
integer walker weights, [sll before recapitulating some 
of the details of the “replica” density-matrix sampling 
technique. [2^, [H, Building upon that previous work, 

our discussion turns to consider the calculation of nuclear 
forces, molecular dipole moments, and atomic dipole po- 
larisabilities, and in so doing confirms the high quality of 
the sampled one- and two-body reduced density matrices 
which is now achievable. 


METHODOLOGY 


i-FCIQMC 


Initiator full configuration interaction quantum Monte 
Carlo provides stochastic integration of the iV-electron, 
imaginary-time Schrddinger equation, yielding wavefunc- 
tions expressed as a linear combination of the set of Slater 
determinants, {|T^i)}, formed from the underlying one- 
particle (most often Hartree-Fock) basis: 


'k = ^Ci|A). (4) 

i 

The coefficients of this wavefunction expansion are ob¬ 
tained by iterative application of the equations 


C, (r + St) = C, {T)-6r (H,, - /r) Q (r)-^ JriJyCj (r) , 

j#i 

( 5 ) 

representing the evolution of the coefficients over a 
timestep St in imaginary time. This evolution is 
achieved by subjecting an ensemble of signed walkers to 
a three-step population dynamics algorithm of “spawn¬ 
ing” , “death”, and “annihilation” steps, the walker pop¬ 
ulations, {IV;}, becoming proportional to the coefficients. 
The full details of this approach have been expounded in 
previous papers, [1S la B 0, EH, E E El , and what 
follows should be regarded only as a brief summary. 

Typically initialised with a single walker placed upon 
the Hartree-Fock determinant, \Dq), a simulation us¬ 
ing integer walkers proceeds with a coupled determinant, 
being randomly selected for each walker on parent 
determinant, |I?i), with a probability Pgen The 

determinant selected, the parent walker then attempts to 
spawn a child on to it, with a probability 


Ps (T>j|A) = 




■^gen 




( 6 ) 


If the attempt is successful, the sign of the spawned 
walker matches that of its parent if iJy < 0 and is in¬ 
verted if Hij > 0. The initiator adaptation, i-FCIQMC, 
modifies the spawning step by introducing a parameter, 
which specifies a lower population threshold under 
which a parent determinant is prevented from spawning 
on to unoccupied determinants. Each walker next at¬ 
tempts to die, with a probability given by 


Pd (A) = St [Ha - , 


( 7 ) 


in which ^ is a population control parameter — known 
as the “shift” — which tends to the ground-state energy 
in the long-r limit. 

These two steps are themselves sufficient to describe 
Eq. [5] fully, but are insufficient to provide convergence to 
a fermionic wavefunction. Instead, a third step — “an¬ 
nihilation” — is required in order to suppress the dele¬ 
terious effects of the fermion sign problem. d. Ell After 
each iteration, walkers of opposite sign on the same de¬ 
terminant annihilate, and in so doing ensure that each 
determinant is populated by walkers of only one sign for 
the next iteration. The success of these processes re¬ 
lies on the sparsity of the wavefunction induced by the 
underlying basis — typically chosen to be Hartree-Eock 
orbitals — which confines it to a generally small region 
of the Hilbert space. In so doing, it ensures that annihi¬ 
lation events are numerous enough to maintain the sign 
structure of the sampled wavefunction accurately. 

Although the walkers of (i)-FCIQMC were initially 
conceived as an ensemble of discrete particles, there 
is some merit in instead positing a set of non-integer 
walkers. [iH 321 Such an approach reduces the amount 


of random number generation required, reduces the in¬ 
stantaneous fluctuations in the populations on a given 
determinant, and hence the fluctuations in the energy 
estimators in imaginary time. 

This formulation is achieved by applying the spawning, 
death, and annihilation steps introduced earlier contin¬ 
uously, rather than discretely. Thus, instead of spawn¬ 
ing a walker of signed integer weight from a determinant 
\Di) to a coupled determinant jUj) with a probability 
Ps (ZIj|iD;), a walker of weight Pg is spawned with prob¬ 
ability 1. Likewise, the death step is remodelled such 
that it simply involves reducing the population on a de¬ 
terminant \Di) by P(j {Di). Annihilation is achieved by 
taking the signed sum of walkers on each determinant on 
a given iteration as the residual population for the next 
iteration. Eor z-FCIQMC calculations, the parameter 
is recast as a continuous variable rather than an integer. 

The continuous nature of the spawned walkers in this 
approach does not, however, imply that the number of 
spawning events becomes continuous. As in the inte¬ 
ger formulation, where there are exactly spawning at¬ 
tempts from determinant \Di) with a population on 
each iteration, there are a discrete number of attempts 
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per determinant per iteration. For practical purposes, 
a continuous spawning threshold, k, is introduced such 
that a Pg < K, K walkers are spawned with a probability 
Ps/K. This implementation is designed to alleviate the 
significant cost of low-weighted spawnings compared to 
their effect on the overall wavefunction, as well as ensur¬ 
ing that the wavefunction remains compact and express¬ 
ible instantaneously by a number of walkers far smaller 
than the size of the space. 

Whilst the death step requires no extra modification 
of this kind, some additional considerations must be ad¬ 
dressed for annihilation. In order that determinants can 
become completely depopulated, and we are not forced to 
store large numbers whose populations are very close to, 
but not exactly, zero, a minimum occupation threshold, 
N^cc^ is imposed upon them. If, after annihilation, the 
population on a determinant its population 

is set either to N^cc with probability Ni/Nacc, or else to 
0 with probability I — Ni/N^cc- 

As a final practical means of alleviating the compu¬ 
tational burden of this approach, it is possible to treat 
only a subspace of the full Hilbert space with non-integer 
walkers, continuing to describe the remainder in a dis- 
cretised fashion. In order to preserve the benefits of the 
non-integer approach on the fluctuations of the energy 
estimators, the truncation is specified by an excitation 
level, X, with only X'fold s-nd lower excitations from the 
reference included in the non-integer subspace. A typical 
choice of parameters N^cc = l)2<x<4(4is used 
here), and k = O.OI entails only a modest increase in the 
computational cost of the calculation over the integer im¬ 
plementation, while retaining many of the benefits of the 
full non-integer approach. 

i-FCIQMC provides two essentially independent en¬ 
ergy estimators, which, taken together, provide a use¬ 
ful confirmation of the validity of the obtained result. 
The first, to which we have already alluded, is the shift, 
/i. This is initially held constant (typically at zero) to 
facilitate an exponential growth in the number of walk¬ 
ers, before being allowed to vary dynamically to keep the 
population constant. At convergence, this variation fluc¬ 
tuates around the energy of system, and thus provides 
an energy on the basis of the growth rate of the entire 
ensemble of walkers. A projected energy estimator, of 
the form 


(i4o|g|d/(T)) 

(^ol'I/(r)) ’ 


( 8 ) 


on the other hand, depends only upon the populations 
of the determinants coupled to the reference state, |I?o). 
Whilst the error in this projection is formally first-order 
in the wavefunction error, its non-variationality tends to 
mean that it converges rather faster to the exact, infinite- 
walker limit than does a variational estimator, owing 
to favourable cancellation of errors. The projected en¬ 
ergy is thus typically preferred when the wavefunction 


is dominated by the Hartree-Fock determinant, but a 
projection on to a multi-reference trial wavefunction or 
the variational estimator provided by the density matri¬ 
ces (which is second-order in the wavefunction error) are 
often more useful in more strongly-correlated cases. [1^ 
Once the ensemble has equilibrated, the simulation is al¬ 
lowed to evolve in imaginary time until the statistical 
errors in both p and Ap^oj have been satisfactorily re¬ 
duced, upon which a Flyvbjerg-Petersen blocking anal¬ 
ysis is performed to estimate the error in the obtained 
result. 371 


Stochastic density-matrix sampling 


In terms of the wavefunction ansatz of i-FCIQMC 
(Eq. |4l) and the creation and annihilation operators, the 
one- and two-body reduced density matrices, 7 and F, 
may be formulated in terms of the wavefunction expan¬ 
sion and the conventional creation and annihilation op¬ 
erators, and a, as 

Ipq = (4'lapSl^) (9) 

= ^CiC'j(Al4«?l^j)> (10) 

ij 


and, 


= (vi,|atat 


apOga^a^l'k) 


/ . (-Hi \ ^p^q^s^r l-Hj ) , 


( 11 ) 

( 12 ) 


respectively, [3^ l39| and an important recent develop¬ 
ment of the theory allows these objects to be sampled 
in an efficient, stochastically unbiased fashion. 2^, 4^ 
The diagonal elements of these objects, of the form 


F 


pqpq 


E 


(13) 


may be calculated straightforwardly, as each determi¬ 
nant, contributes Cf to each of the matrix 

elements involving its occupied orbitals. The correspond¬ 
ing explicit generation of all the required determinant 
pairs for the off-diagonal elements is not practical, but 
the observation that the relevant pairs are at most dou¬ 
ble excitations of one another allows both 7 and F to be 
sampled via the spawning steps, Thus, the existing 
computational effort required for the communication of 
the spawning event need only be slightly accentuated (by 
the need now to convey both the amplitude and the iden¬ 
tity of the parent determinant to the child) to allow the 
contributions to the off-diagonal matrix elements from 
determinant pairs to be calculated on the fly. 

As these off-diagonal contributions are only added 
upon a successful spawning event, it is necessary that 
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they be rescaled according to the probability of such an 
event taking place. That is, a contribution CjCj will in¬ 
stead be accumulated as 


CiCj ( D, 


apa^aga^ 




Pc (-DjIA) 


(14) 


with (A|.C);) the probability that at least one spawn¬ 
ing attempt from \Di) to |A) is successful on a given 
iteration. Depending on whether integer or non-integer 
walkers are considered, this probability is given by 


1 - 


A e z 

n,^z, 

(15) 


with A fhs instantaneous walker population residing on 
iDj) and A the probability that no walker is spawned 
between \Di) and | A) in a single attempt. For an integer 
spawning event, this probability is 

Aint = 1 - min {6t |Aj| ,Pgen (AIA)) , (16) 

but this must be modified in the case of continuous 
spawning to 


in terms of the one- and two-electron integrals, {Ag} 
and I Qpqrs }; as 


^pqrs ‘2N 2 ^qs^pr) 9pqrs: ( 1 ^) 

will give rise to an biasing error in density matrices in the 
long-T limit for determinant pairs where fcpgr-s ~ 6, but 
whose amplitudes are both significantly non-zero. How¬ 
ever, modifications to the algorithm to treat the bias re¬ 
maining beyond that already defined by Brillouin’s theo¬ 
rem explicitly — such as introducing additional events to 
spawn walkers proportionally to the inverse of the hamil- 
tonian element — have been shown to be of little addi¬ 
tional benefit due to the negligible nature of this bias in 
numerical studies to date.[29| 

In a separate difficulty, it has been shown previously 
that a straightforward implementation of the above sam¬ 
pling gave rise to a convergence of the density matrices 
with increasing Av which was rather slower than that 
of, say, the projected energy. This behaviour stems not 
simply from undersampling, but rather from a bias in 
the statistical sampling technique itself. In particular, 
appropriate contributions to the matrix elements are ap¬ 
proximated by 


A 


cont 


1 - 




.1 - Pgen (A I A) 


Pg< K 

otherwise, 


(17) 


where k is the continuous spawning threshold, if used. 

A naive implementation of the above sampling is satis¬ 
factory for the accumulation of approximate density ma¬ 
trices, but is beset by a number of shortcomings which 
should be considered. [29j| As contributions to the off- 
diagonal matrix elements are only added upon a suc¬ 
cessful spawning attempt, problems can arise when the 
spawning events are discretised. In this case, the proba¬ 
bility that such an event occurs is proportional to the cou¬ 
pling hamiltonian matrix element, Aj: pairs of de¬ 

terminants which are connected by large matrix elements 
are correspondingly sampled more often than pairs which 
are only weakly coupled. Thus, if two highly-weighted 
determinants contained in the stochastically-sampled, in¬ 
teger walker space are connected by a small hamiltonian 
element, their contribution to the density matrices may 
be severely under-represented, or neglected entirely. 

This problem is most notably in evidence in the case of 
single excitations of the reference determinant, for which 
the coupling matrix elements are strictly zero accord¬ 
ing the Brillouin’s theorem. This is countered in the 
present implementation by accounting for these contribu¬ 
tions to the density matrices explicitly, and hence remov¬ 
ing the dependence upon a successful spawning event. 
Other contributions, however, whose sampling will still 
be proportional to the reduced hamiltonian. |4l| defined 


(A (r)).(A (t)). = (A(r)A(^)). - a(A(T), A(^)) 

(19) 

^(A(r)A(^)>.> (20) 


ignoring the potentially significant covariance, a, be¬ 
tween the two amplitudes and introducing a bias, 
whether or not the averaged walker populations are 
themselves unbiased. It is, to that end, unsurprising that 
this problem is at its greatest for diagonal elements, for 
which the “two” amplitudes are perfectly correlated, and 
— the error being of a single sign — there is no possibility 
of error cancellation. 


This problem is rather more serious than the previous 
concerns over discretised spawning, but one for which a 
rather simple solution exists. Unbiased density matri¬ 
ces can be calculated with the introduction of a second, 
uncorrelated walker ensemble, to which the stochastic 
spawning, death, and annihilation steps are applied inde¬ 
pendently, and whose statistics are acquired separately, 
from the first.This adaptation, known as replica sam¬ 
pling, achieves the unbiasing by ensuring that all the 
products of determinant amplitudes are calculated us¬ 
ing populations from both simulations, and has previ¬ 
ously found application in the stochastic sampling of 
the A^-electron density matrix known as density matrix 
quantum Monte Carlo, 42j and the recently-introduced 
Krylov-projected quantum Monte Carlo, That is, for 
example, a successful spawning event from \Di) to |A) 
in replica 1, occurring with a probability (A|.C);), 
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gives rise to a contribution of: 


^a)^,(2) 


^(2)^(1) 


+ 


,(2) 


(^jIA) 


( 21 ) 


This approach bears some conceptual similarity with the 
bilinear sampling algorithm in Green’s function Monte 
Carlo, introduced by Zhang and Kalos, in that both seek 
a means of finding expectation values of operators which 
do not commute with the hamiltonian, via two sets of 
independent walker distributions, [s^ The main differ¬ 
ence, though, is that the bilinear approach transforms 
the Schrodinger equation such that there are two related 
wavefunctions to sample, while in the present work the 
walker ensembles are independent samples of the same 
underlying object. In providing a stochastically unbi¬ 
ased route to the density matrices, the replica sampling 
technique thus provides the first realistic opportunity to 
achieve high-accuracy ab initio results for the sizeable 
suite of properties that can be derived therefrom. 


NUCLEAR FORCES 


The force acting on a nucleus in a molecule or cluster is 
defined as the negative gradient of the molecular energy 
with respect to the nuclear coordinates: 


F = - 


dE 


( 22 ) 


In Eq. (1^ the symbol F denotes the nuclear force vec¬ 
tor, E the energy of a molecule at a fixed geometry in 
the electronic ground state, and R refers to the vector of 
nuclear coordinates in the centre of mass frame of refer¬ 
ence. A comprehensive review of techniques and explicit 
expressions to compute derivatives of the electronic en¬ 
ergy with respect to nuclear coordinates is available in the 


literature. [4,11 - 1451 The following discussion is thus limited 


to the basic concepts for the calculation of nuclear forces 
using all-electron FCI wavefunctions as obtained as a sta¬ 
tistical average using the z-FCIQMC method, once the 
calculation has been converged with respect to the num¬ 
ber of walkers. In the present work, we have adjusted 
the total number of walkers to achieve a walker popu¬ 
lation of 50000 at the reference (i. e. highest populated) 
determinant. Preceding work confirmed that, at such 
a population levels, noise arising from small stochastic 
populations of random determinants is sufficiently sup¬ 
pressed and the wavefunction converged. 

The first derivative of the electronic energy of a Cl 
wavefunction generally depends on the derivatives of the 
atomic orbitals (AOs) and the molecular-orbital (MO) 
and Cl coefficients. All these terms depend upon the nu¬ 
clear coordinates, and the computation of nuclear forces 
requires knowledge of the first derivatives with respect 
to all considered degrees of freedom. However, electronic 


wave functions obtained from i-FCIQMC optimizations 
are variational with respect to the Cl coefficients and 
a component F^. of the nuclear force vector can be ex¬ 
pressed in terms of the reduced density matrices as 


^ dhpq r d{pq\rs) 

dx dx 

pq pqrs 


in which the terms {hp^} represent the one-electron in¬ 
tegrals from the hamiltonian, and is the contribution 
from the fixed nuclei. Moreover, all-electron FCI wave- 
functions considered in this work are also invariant un¬ 
der variation of the MO coefficients. The nuclear forces 
can thus be expressed solely in terms of the one- and 
two-electron density matrices and the skeleton derivative 
integrals of the basis functions: 


MO AO Q, 

F =-VV7 c c 


pq py 

MO AO 


~ ^ ^ ^pqrs^^ip^^vqCprC, 


dx dx 

d{pv\p(T) 


pqrs pyp<7 

MO AO 


dx 




dS„ 


py 

d~’ 


(24) 


where 


MO 


MO 


Xpq E ^pr^qr ^ E ^prst {qr\st) (25) 


is an element of the lagrangian and is an element of 
the overlap matrix. In particular, neither the computa¬ 
tion of derivatives of the Cl hamiltonian matrix nor the 
solution of the coupled-perturbed Hartree-Fock equation 
for the derivatives of the MO coefficients are required. 
We have implemented an interface to MOLPRO to compute 
the integrals and nuclear forces from the f-FCIQMC den¬ 
sity matrices, [i^ 

As a first benchmark, we have applied the z-FCIQMC 
methodology to compute the nuclear forces at several 
points along the dissociation curve of molecular nitro¬ 
gen, as the electronic wavefunction changes from single- 
to strong multi-reference character. Figure [T] (top) com¬ 
pares the potential energy computed with i-FCIQMC 
and the FCI program in MOLPRO using a small 6-3 IG ba¬ 
sis set to allow for comparison to exact (FCI) results. 
The accuracy of the i-FCIQMC methodology for the 
computation of total energies was already evaluated 0, 
and we generally find excellent agreement between the 
i-FCIQMC and FCI data set. 

In Figure [T] (bottom), the nuclear forces for the same 
geometries are illustrated. Comparison with FCI results 
obtained from numerical gradients provides a direct mea¬ 
sure of the quality of the reduced density matrices com¬ 
puted from the replica algorithm based on f-FCIQMC, 
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FIG. 1. Top: Potential energy profile for the N-N bond dis¬ 
sociation of N 2 relative to the energy of two isolated nitrogen 
atoms in the electronic ground state. Bottom: corresponding 
forces at one nitrogen atom computed using analytic gradients 
from i-FCIQMC reduced density matrices, compared to FCI 
with numerical differentiation. Results are identical within 
the accuracy of the numerical differentiation. The respective 
minimum energy = -0.2685 a.u.) and force (F^jin = 

0.0 a.u.) at an internuclear distance of 2.144 a.u. is indicated 
by the blue symbols. All results were obtained with a 6-31G 
basis set. 


and, once again, the data shows excellent agreement be¬ 
tween the analytic i-FCIQMC forces and the FCI results 
for all geometries. 

As second example for the calculation of analytic gra¬ 
dients and nuclear forces, we considered symmetric dis¬ 
placements of the atoms in a water molecule along the 
OH bonds. In a small 6-31G basis set, exact (FCI) di- 
agonalisation of the hamiltonian matrix is still feasible 
and Figure [2] illustrates results from FCI reference and 
i-FCIQMC calculations. The nuclear forces as shown in 
Figure [2] have been obtained from the Cartesian force 
vectors as the absolute force acting on either a hydro¬ 
gen or the oxygen atom with the sign taken from the z- 
component of the force vector, which has been aligned 
with the symmetry principal axis. Although there is 
no computational advantage over direct diagonalisation 
methods for basis sets as small as the 6-31G basis, the 
replica algorithm implemented in i-FCIQMC can be ap¬ 
plied to much larger molecules and basis sets, providing 
essentially numerically exact nuclear forces. In order to 
demonstrate the scope of the i-FCIQMC replica technol¬ 
ogy, we have also computed the all-electron forces within 
a cc-pVTZ basis set, evidently an infeasible task for cur¬ 
rent deterministic FCI algorithms, where the many-body 
basis now spans (!I[10^^] determinants. Figure [5] (dashed 
lines) illustrates the notably larger forces at intermediate 


FIG. 2. Absolute forces acting on the oxygen and hydrogen 
atoms in a H 2 O molecule computed using i-FCIQMG and FCI 
with a 6-31G and cc-pVTZ basis set (the sign corresponds to 
the z-component of the force vector). The data were acquired 
for symmetric displacements of the hydrogen atoms from the 
equilibrium geometry. The abscissa indicates the OH bond 
length of the respective molecular geometry. 

stretching of the OH bonds if accurate cc-pVTZ basis set 
are combined with this level of theory in the calcula¬ 
tions. This would have implications for dynamics calcu¬ 
lations, as well as providing the basis for highly accu¬ 
rate geometry optimisations for systems with electronic 
ground states of strong multi-reference character. 


THE DIPOLE MOMENT OF CO 

The interaction of an electronic system of charge q with 
an external electric field, in an external potential, V, 
may be expressed as an expansion in terms of multipoles, 

E = qV-^.^-le-^-..., (26) 

with /X the rank-1 dipole moment, © the rank-2 
quadrupole moment, and so on. It is the dipole moment 
itself with which we are presently concerned, and which 
may be calculated according to: 

M = (d/1 A I d/) (27) 

N 

= (d/|^9,r,|vl/) (28) 

i 

N 

= -^(vI/|rJvl,), (29) 
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where, in the last line, the substitution g, = — 1 (for 
electrons) has been made. Applying the Slater-Condon 
rules, [13, 13 this expression can be recast in terms of 
the one-body reduced density matrix and one-electron 
molecular-orbital integrals for an arbitrary Cartesian 
component, w, as, 

= -^lpq{(l>p\w\(j)q) +'^ZjR^^\ (30) 

pq I 


to which the contribution from the (fixed) nuclei with 
charges {Z^} and positions {Rj} has been added. Thus, 
given the molecular-orbital integrals, x \ (pq)}, 

|(d>r, 11/1 and {((pp \ z \ (pq)}, which are readily 

available, [46l 13l the one-body reduced density matrix 
obtained from i-FCIQMC provides direct access to the 
dipole moment, and, more generally, to multipole mo¬ 
ments of arbitrary rank. 

As an interesting application of this approach, we con¬ 
sider the well-known problem of the dipole moment of 
CO at its equilibrium bond length, 2.1316 OQ.jsOl This 
system, with its subtle combination of a and tt effects, is 
difficult to predict intuitively a priori, and Hartree-Fock 
theory notably suggests the polarity to be C''’0~, while 
it is experimentally known to be C^O'*'. 

We use the large aug-cc-pVATZ-DK basis sets for 
this study and adopt the second-order Douglas-Kroll- 
Hess hamiltonian. [SlMSSj Although relativistic effects are 
small for these comparatively light atoms, the calculation 
of the dipole moment tends to be strongly basis-set de¬ 
pendent, and the use of a large set becomes correspond¬ 
ingly desirable. To that same end, it is desirable to be 
able to extrapolate finite-basis dipole moments to the 
complete-basis-set limit, as such extrapolations have pre¬ 
viously been useful in AFCIQMC studies. 0 It has been 
shown that the asymptotic convergence of the correla¬ 
tion part of the dipole moment with the cardinality of 
the basis set, X, is suitably described by the form 


„(X) _ (CBS) 

H'corr H'corr 


aX 


-3 


(31) 


in much the same way as the correlation energy itself, 

57l| The complete-basis-set limit correlation contribution 
to the dipole moment, may thus be derived from 

two consecutive finite-basis results, of cardinality X — 1 
and X, according to 


.1^3 (X) 1^3 (X-1) 

(CBS) _ ^ Mcorr -^1 Mcorr 

- (X - 1)3 


(32) 


to which the Hartree-Fock contribution in a suitably 
large basis (aug-cc-pV5Z-DK is used here, for which 
Mz.hf = —0.10355 eoo) may then be added to obtain the 
total dipole moment. 

Table H] presents the results of this approach, with the 
extrapolations performed from the double- and triple-i^ 
and the triple- and quadruple-^ basis sets, alongside the 


Mz/eao 



aug-cc-pVXZ-DK 

CBS 


X =D 

X =T 

X =Q 

(DT) 

(TQ) 

HF 

-0.10135 

-0.10435 

-0.10369 

- 

- 

MRCI 

0.07175 

0.07203 

0.07066 

0.07419 

0.06929 

CCSD 

0.06829 

0.05594 

0.05087 

0.05278 

0.04681 

i-FCIQMC 0.05893(3) 0.05200(4) 0.0474(4) 0.05112 

0.0437 


TABLE I. Calculated dipole moments, for CO at the 
HE, MRCI (using a 10-electron, 8-orbital active space), [s^- 
IIO] CCSD.[Sl] and i-FCIQMC levels of theory, with the 
complete-basis-set limit obtained from two-point, inverse- 
cube extrapolations. [13] The standard error (in brackets) is 
derived as the standard deviation of the results from three in¬ 
dependent i-FCIQMC calculations. The experimentally ob¬ 
tained bond length, 2.1316 Oq, is used.fs^j and the la la*^ 
electrons are held frozen and neither relaxed nor optimised 
for the response of an electric field. The signs are arranged 
such that < 0 indicates a C^O~ polarity, and thus all 
the post-Hartree-Fock methods successfully reproduce qual¬ 
itative agreement with the observed dipole’s direction. The 
i-FCIQMC calculations were performed for 24 hours on 400 
cores {X =D and T) or 600 cores (A =Q) using 


walkers, with the adjustable parameters = 1, x = 4, 

K = 0.01, and = 3.0, and the timstep allowed to vary 
dynamically to limit noisy walker growth. The sizes of the 
full orbital spaces for the double-, triple-, and quadruple-^ 
calculations are 44, 90, and 158, respectively. 


analogous results from coupled-cluster theory and multi¬ 
reference CL The rapid convergence of with number of 
walkers in the i-FCIQMC dynamic is also demonstrated 
in Figure [31 



FIG. 3. Calculated dipole moments for CO in an aug-cc- 
pVDZ-DK basis as a function of the number of walkers. In¬ 
creasing the walker population is beneficial in reducing the 
stochastic error in the final result, but a qualitative descrip¬ 
tion of the system is achieved at rather modest N^, as indi¬ 
cated by the fineness of the scale presented here. 


By comparison with the experimental dipole moment, 
variously given as 0.044 eug and 0.048 ean. |6ll464l| it is 












{ji^leaQ Abs. relative error (%) 


apparent that i-FCIQMC performs rather better than 
MRCI, and is comparable to CCSD. However, it can be 
seen that CCSD actually overestimates the dipole mo¬ 
ment compared to i-FCIQMC, which can be taken as 
close to exact in each of the finite basis sets, and this 
feature of CCSD allows for favourable cancellation of er¬ 
rors with the basis-set incompleteness, yielding the fortu¬ 
itously accurate extrapolated result. The remaining dis¬ 
parity between these results and experiment should not 
be ascribed to an inadequacy of the z-FCIQMC density 
matrices, but is rather largely attributable to basis-set 
incompleteness error. Indeed, the larger (TQ) extrapo¬ 
lations are rather more satisfactory than the correspond¬ 
ing (DT) results, highlighting the sensitivity of such ap¬ 
proaches to the adequacy of the choice of basis. This 
effect has been previously observed in the context of ion¬ 
isation potentials, 12| but is magnified in this instance 
by the stronger basis-set dependence of dipole moments 
than correlation energies. It is also worth noting that 
a small vibrational contribution to the dipole moment 
is expected.but the results of this study support the 
view expounded in that work by Luis and coworkers, that 
an accurate treatment of electron correlation in a suffi¬ 
ciently large basis set is adequate for close agreement 
with experiment. 


The quality of the i-FCIQMC density matrices may be 
illustrated by considering the CO problem in a small cc- 
pVDZ basis, for which deterministic FCI results can be 
obtained. In this case, whose results are summarised in 
Table mi i-FCIQMC reproduces the FCI dipole moment 
to within 0.06%, whilst the CCSD and MRCI results are 
in error by 10% and 8%, respectively. Also of note is 
that, whilst the quoted i-FCIQMC result was obtained 
using (^0®^ walkers, it can be obtained just as well, 
and with apparently negligible initiator error, with only 

These results, therefore, bear out the supposed 
high quality of the sampled density matrices, and in 
demonstrating the compatibility of i-FCIQMC with the 
Hellmann-Feynman theorem, suggest that future studies 
of energy derivatives and their associated properties may 
well prove fruitful. 


ATOMIC DIPOLE POLARISABILITIES 


The previous section began by noting the dependence 
upon the permanent dipole moment of a system’s inter¬ 
action with an applied electric field as given by — 

Of course, the application of such a field will, in reality, 
affect the distribution of charge, and hence the dipole mo¬ 
ment itself. Expanding the dipole moment as a function 
of the field, therefore, we may write a given component. 


HF 

-0.0915 

201.10 

MRCI 

0.0973 

7.61 

CCSD 

0.0996 

9.94 

CCSDT 

0.0931 

2.87 

CCSDTQ 

0.0906 

0.11 

CCSDTQP 

0.0905 

- 

i-FCIQMC 0.09045(3) 

0.06 

FCI 

0.0905 

- 


TABLE II. Comparison of obtained dipole moments of CO in 
a small cc-pVDZ basis to the deterministic FCI result. As in 
Table m at all levels of theory, the two core orbitals were held 
frozen, and neither relaxed nor optimised for the response of 
an electric field. The i-FCIQMC result being in error by less 
than 0.1%, the density matrices derived therefrom are thus 
shown to be of near-FCl quality. The coupled-cluster results 
— obtained by finite differentiation (±2 x 
using the MOLPRoji^, [o^ and MRCc[63] codes — are slow to 
converge to the FCI limit, with quadruple excitations needed 
for high accuracy. 


/ij, as 

Ml = ^ + 2 ^ + • ■ •, (33) 

3 jk 

where a^j and Pijj. represent elements of the polarisabil- 
ity and first hyperpolarisability tensors, respectively. 

is the zero-field permanent dipole, which is always 
zero for atomic species. Whilst the effect of the induced 
dipole moment is generally less significant for polar sys¬ 
tems, it is the leading-order term in the expansion of 
the dipole moment for atoms which has no static dipole. 
It is thus crucial in accounting for the dipole-dipole dis¬ 
persion interactions which often bind such species, and 
indeed will be the first-order response not only to static, 
but also to dynamic fields, The calculation of ol thus 
provides an interesting study in and of itself, as well as 
a probing test of the calculation of reduced density ma¬ 
trices with i-FCIQMC. We here consider the noble-gas 
atoms, Ne, Ar, and Kr, as archetypal examples of the 
problem. 

It is apparent from Eq. [33] that the polarisability may 
be thought of as the derivative of the dipole with respect 
to the field. 




ry ■ ■ =. 


«=o 


(34) 


evaluated at ^ = 0. As for many response properties, 
this may be calculated by solution of the coupled per¬ 
turbed Hartree-Fock equations, [ t^ but for our purposes 
it is convenient to suppose that a particular component, 
Qfjj, might be effectively calculated by a finite-difference 
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/ 2 2 r ^-1 

a^z/e flo-Eh 
aug-cc-pVTZ-DK 


System 

Analytic 

Finite-field Abs 

1 . relative error 

Ne 

2.438384 

2.437962 

1.73xl0"'‘ 

Ar 

10.841398 

10.842952 

1.43xl0"'‘ 

Kr 

16.674939 

16.680012 

3.04x10""^ 



aug-cc-pVQZ-DK 


Analytic 

Finite-field Abs 

1 . relative error 

Ne 

2.620174 

2.619806 

1.40xl0”'‘ 

Ar 

11.128300 

11.131348 

2.74x10""^ 

Kr 

16.792296 

16.799500 

4.29x10""^ 


TABLE III. Analytic MP2 dipole polarisabilities, for the 
noble gases Ne, Ar, and Kr, in aug-cc-pVTZ-DK and aug-cc- 
pVQZ-DK basis sets compared with the corresponding finite- 
field results, calculated with an electric field strength of 0.005 
-Eh/eap. 


approach, 


At, (<5^,) -,, , 

in which is a small field applied in the j direction, 
and Hi (S^j) is the ith component of the dipole moment 
induced by so doing. The second equality holds for the 
spherically-symmetric atomic systems under considera¬ 
tion here since = 0 , and the errors resulting formula 
are second-order since it is now equivalent to a central- 
difference approximation. 

Straightforward and appealing though this implemen¬ 
tation is, it is useful before proceeding to have some no¬ 
tion of its performance relative to analytic gradient meth¬ 
ods. In particular, analytic gradients are readily and 
rapidly available from MP2 theory, and this approach 
thus provides a useful framework in which to assess the 
suitability of the finite-held method. 

The results of this comparison, with the finite-field 
polarisabilities performed in a field of strength 0.005 
E^^jeaQ, are summarised in Table IIIII The mean abso¬ 
lute percentage error inherent in the approach is found 
to be of the order of 0.02%, demonstrating both its suit¬ 
ability for the purpose, and also that the field chosen 
is sufficiently small to establish the pseudo-linear depen¬ 
dence of the induced dipole upon the field. That this 
dependence is established without having to use a very 
small field is encouraging, since in the stochastic formu¬ 
lation provided by f-FCIQMC, the stochastic error in the 
induced dipole must be divided by the field strength to 
obtain the equivalent error bounds in the polarisability. 
This behaviour is illustrated for Ne in the aug-cc-VTZ- 
DK basis in FigureSl which highlights the balance which 
must be achieved between minimising second-order ef¬ 


fects and maintaining a suitable level of stochastic error. 



FIG. 4. Calculated dipole polarisabilities for the Ne atom 
in an aug-cc-pVTZ-DK basis with different applied field 
strengths, As in the previous section, the i-FCIQMC cal¬ 
culations were performed using G (^10®^ walkers, a dynamic 
timestep, and the adjustable parameters = 1, y = 4, 
K = 0.01, and = 3.0. Sufficiently small fields establish the 
required pseudo-linear relationship between the polarisability 
and the applied field, but too small a field gives rise to large 
stochastic errors. This initial study prioritises the elimination 
of non-linear effects, and a field strength of = 0.005 E^^jea^ 
is thus chosen as suitable for the remainder of this work, where 
the random errors can be systematically controlled. It is en¬ 
couraging to note, however, that this choice may be somewhat 
conservative, and that a slightly larger field may be permissi¬ 
ble in future work. 


Secure in the knowledge of the suitability of the finite- 
field approach, we may proceed with an assessment of the 
performance of f-FCIQMC compared with other meth¬ 
ods. Specifically, as in the previous section, we cal¬ 
culate the polarisabilities using CCSD and MRCI for 
comparison, though in this case the extrap¬ 

olations are performed from results at the triple- and 
quadruple-C basis sizes, reflecting the increased sensitiv¬ 
ity to basis set incompleteness error which this quantity 
entails. 


As might have been expected, the error incurred by 
extrapolating is somewhat reduced upon including the 
larger quadruple-^ treatment, and the i-FCIQMC results 
given in Table IIVI bear correspondingly close agreement 
with experiment. The remaining errors — in the region 
of 0.5 to 3% — are nonetheless still likely to be artefacts 
of the basis sets, as the application of a field accentu¬ 
ates the importance of describing the intricacies of the 
more diffuse regions of electron density. Thus, although 
“augmented” basis sets are employed, there is likely still 
something to be gained from a more complete description 
of this behaviour. This suggestion is, once again, further 
strengthened by the fact that i-FCIQMC is able to re¬ 
cover FCI-quality results for basis sets in which direct 
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/ 2 2 r ^-1 

System aug-cc-pVTZ-DK aug-cc-pVQZ-DK CBS Experiment 


Ne 

2.42(1) 

2.596(2) 

2.65 

2.57 

Ar 

10.855(5) 

11.092(3) 

11.08 

11.23 

Kr 

16.81(4) 

16.86(6) 

16.82 

16.73 


TABLE IV. i-FCIQMC polarisabilities of the noble gases Ne, 
Ar, and Kr, obtained in aug-cc-pVTZ-DK and aug-cc-pVQZ- 
DK basis sets, along with the extrapolated complete-basis-set 
limit results. The number in brackets indicates the error in 
the preceding digit, obtained as the standard deviation of the 
results of three independent calculations. The experimen¬ 
tal results are shown for comparison. [tH, The i-FCIQMC 

calculations were performed using G (^10®^ walkers, the ad¬ 
justable parameters = 1, x = 4, k = 0.01, and = 3.0, 
and run for 48 hours on 320 cores. 

/ 2 2 n-l 

a^^/e floEh 

aug-cc-pVTZ-DK aug-cc-pVQZ-DK CBS 


System 

HE 

CCSD 

MRCI 

HE 

CCSD MRCI CCSD MRCI 

Ne 

2.20 

2.42 

2.41 

2.33 

2.59 

2.58 

2.64 

2.64 

Ar 

10.45 

10.81 

10.36 

10.72 

11.03 

11.53 

11.00 

12.20 

Kr 

16.21 

16.78 

17.15 

16.36 

16.81 

17.22 

16.74 

17.18 


TABLE V. Polarisabilities of the noble gases Ne, Ar, and Kr, 
computed using Hartree-Fock, coupled-cluster, and multi¬ 
reference Cl (with an 8-electron, 8-orbital active space) the¬ 
ories. As in Table the extrapolations to the complete- 
basis-set limits are also shown. [tD. It^] 


comparison is possible, reproducing the polarisability of 
Ne in a small cc-pVDZ basis to within 0.005%, for in¬ 
stance. 

The same results, computed using Hartree-Fock the¬ 
ory, CCSD,[i^ and MRCI, are listed in Table Ivl 

The mean (absolute) error for the MRCI calculations is 
4.7%, whilst that for i-FCIQMC, and coupled-cluster the¬ 
ory, is around 1.6%. The comparability is unsurprising, 
given the ascription of much of the error to finite-basis 
effects. However, it is now necessary to investigate the 
impact of stronger correlation on this quantity in more 
challenging systems, where we expect more significant 
advantages to come from i-FCIQMC. 


CONCLUSIONS 

The results presented in this work serve to confirm 
the high quality of the stochastically-obtained reduced 
density matrices available via replica sampling in i- 
FCIQMC, capable as they are of reproducing FCI-quality 
results for nuclear forces, dipole moments, and polaris¬ 
abilities, and in some cases close agreement with exper¬ 
imental values. In so doing, they cement the place of 
the replica technique as an important extension to the 


theory, and widen its scope considerably. 

In addition to the most obvious extension of an abil¬ 
ity to compute a larger range of properties for a wider 
variety of systems, there remain a number of theoretical 
and technical challenges to be addressed in future stud¬ 
ies. Perhaps the most pressing task is to extend this 
work to encompass results from open-shell systems, in 
which correlation effects are likely to be more important. 
Moreover, if comparisons to experimental results are to 
be further sought and achieved for dipole moment prop¬ 
erties, there is some motivation to explore larger basis 
sets with multiple levels of augmentation. [731 l74l| which 
may be of particular use in better describing the more 
diffuse electron densities of finite-field calculations, and 
more generally in describing larger and heavier atoms of 
interest. 
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